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ABSTRACT 


This  survey  of  image  restoration  techniques  presents  a  concise  overview  of 
the  most  useful  restoration  methods.  Linear  spatially  invariant  and  linear  spa¬ 
tially  variant  image  restoration  techniques  are  described  and  the  strengths  and 
weaknesses  of  each  approach  are  identified.  Examples  of  restored  images  for 
the  various  techniques  are  given.  To  provide  guidelines  for  choosing  a  restora¬ 
tion  technique  for  a  particular  application,  a  comparison  of  the  techniques  is 
made.  The  restoration  methods  are  compared  and  evaluated  based  on  the  fol¬ 
lowing  criteria:  restored  image  visual  quality,  performance  in  the  presence  of 
additive  image-independent  noise,  degree  of  a  priori  knowledge  required,  and 
computational  complexity. 
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EXECUTIVE  SUMMARY 


This  technical  report  gives  a  broad  overview  of  the  field 
of  image  restoration  and  focuses  on  the  most  commonly 
used  and  successful  restoration  techniques.  Comparisons 
of  these  techniques  are  made  so  that  the  reader  who  is 
faced  with  a  specific  application  can  quickly  perform  the 
necessary  trade-offs  to  choose  an  appropriate  restoration 
method. 

Image  restoration  is  defined  as  the  processing  per¬ 
formed  on  a  degraded  image  to  remove  or  reduce  the  de¬ 
grading  effects.  Restoration  differs  from  image  enhance¬ 
ment  in  that  it  requires  prior  knowledge  of  the  degrada¬ 
tion  phenomena.  The  degradations  in  an  imaging  system 
arise  from  the  filtering  effects  of  the  optical  elements, 
electrical  components,  and  surrounding  environment. 

To  lay  the  foundation  for  the  analytic  methods  dis¬ 
cussed,  basic  models  of  imaging  systems  are  introduced 
in  Section  2.0  and  the  effects  of  these  systems  on  the 
original  object  are  identified.  Linear  spatially  variant  and 
invariant  as  well  as  spatially  separable  systems  arc  de¬ 
fined.  Since  restoration  techniques  require  knowledge  of 
the  degrading  phenomena,  methods  for  determining  these 
degradations  a  posteriori  from  an  output  blurred  image 


are  given  in  Section  3.0.  Methods  of  direct  measurement 
of  the  system  point  spread  function  from  the  images  of 
point,  line,  and  edge  sources  are  outlined.  Estimation 
techniques  that  use  homomorphic  signal  processing  con¬ 
cepts  are  also  presented. 

The  main  body  of  the  report,  Section  4.0,  examines 
the  derivation  of  linear  spatially  invariant  and  variant 
restoration  techniques.  The  restoration  criterion  motivat¬ 
ing  each  technique  is  supplied  and  the  assumptions  and 
limitations  of  each  method  are  identified.  Techniques  de¬ 
rived  include  simple  inverse  filtering,  constrained  decon¬ 
volution,  homomorphic  deconvolution,  and  maximum 
a  posteriori  estimation.  Superresolution  methods  are  also 
briefly  discussed. 

Section  5.0  compares  the  image  restoration  techniques 
considered.  The  methods  are  compared  against  the  fol¬ 
lowing  criteria:  visual  restoration  quality,  noise  perfor¬ 
mance,  a  priori  knowledge,  and  computational  complex¬ 
ity.  The  information  in  this  section  offers  guidelines  for 
choosing  a  restoration  technique  for  a  particular  applica¬ 
tion.  Conclusions  and  potential  areas  for  future  investiga¬ 
tion  are  discussed  in  Section  6.0. 
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1.0  INTRODUCTION 


This  technical  memorandum  surveys  the  field  of  im¬ 
age  restoration  and  lays  the  theoretical  groundwork  need¬ 
ed  to  understand  the  concepts  and  alternatives  in  this 
branch  of  image  processing.  Although  the  scope  of  the 
report  is  general  in  nature,  the  reader  with  a  specific  ap¬ 
plication  should  be  able  to  quickly  compare  and  contrast 
the  restoration  techniques  applicable  to  his  or  her 
problem. 

1.1  DEFINITION  OF  RESTORATION 

Image  restoration  is  defined  as  the  processing  per¬ 
formed  on  the  image  using  a  priori  knowledge  of  the 
degrading  phenomena  to  reduce  or  remove  the  degrad¬ 
ing  effects.  Restoration  differs  from  image  enhancement, 
which  uses  primarily  ad  hoc  techniques  (frequently  non¬ 
linear)  to  bring  out  or  alter  certain  aspects  of  the  image. 
Enhancement  techniques  include  histogram  reshaping  for 
gray-level  modification,  “crispening"  to  sharpen  edges, 
and  thresholding  to  segment  the  image.  Image  restora¬ 
tion,  on  the  other  hand,  seeks  to  undo,  partially  or  whol¬ 
ly.  the  degradation  an  image  has  undergone  for  the  pur¬ 
pose  of  restoring  it  to  its  original  ideal  form. 

To  undo  the  degradation,  we  must  have  prior  knowl¬ 
edge  of  the  degrading  process.  The  a  priori  information 
required  by  a  restoration  technique  may  simply  be  a  gener¬ 
al  analytic  model  for  the  degrading  imaging  system.  An¬ 
other  method  may  require  detailed  information  about  the 
system  filtering  effects  and  the  imaged  object  and  noise 
characteristics.  The  degree  of  prior  knowledge  available 
for  a  specific  application  will  thus  influence  the  choice 
of  restoration  technique.  In  Section  4  0  we  consider  resto¬ 
ration  techniques  that  vary  widely  in  the  degree  of  prior 
knowledge  required. 

Major  sources  of  image  degradation  and  specific  resto¬ 
ration  examples  are  identified  in  this  introductory  section. 
The  general  analytic  imaging  system  models  used  to  derive 
the  restoration  approaches  are  presented  in  Section  2.0. 
Section  3.0  discusses  the  a  posteriori  determination  of  the 
imaging  degradations  by  direct  measurement  and  estima¬ 
tion  methods.  This  information  is  subsequently  used  as 
the  prior  knowledge  requisite  for  the  restoration  tech¬ 
niques.  Both  linear  spatially  invariant  and  variant  tech¬ 
niques  are  outlined  in  Section  4.0,  along  with  a  few  non¬ 
linear  restoration  methods.  Section  5.0  offers  a  compari¬ 
son  of  the  techniques  (Table  3  is  especially  useful)  and 
criteria  to  consider  when  choosing  a  restoration  method 
for  a  particular  application.  Conclusions  and  additional 
comments  arc  given  in  Section  6.0. 


1.2  SOURCES  OF  IMAGE  DEGRADATION 

In  a  real  imaging  system,  the  image  is  degraded  by  the 
net  effect  of  numerous  degradation  phenomena.  Individu¬ 
al  degradations  (other  than  system  noise)  are  typically 
modeled  as  having  a  linear  filtering  effect  on  the  original 
object.  In  general,  the  individual  degradation  effects  may 
interact  in  a  complicated  fashion  to  produce  the  net  effect 
on  the  image.  However,  the  overall  system  degradation 
is  usually  modeled  as  the  linearly  cascaded  effects  of  the 
individual  degradations  present  in  the  system  to  ensure 
analytic  tractability. 

Image  degradations  arise  from  a  variety  of  sources,  such 
as  imaging  system  optics  and  electronics,  and  surrounding 
environmental  effects.  Several  major  sources  of  image  de¬ 
gradation  are: 

1 .  Image  motion.  When  the  object  or  imaging  system 
experiences  relative  motion  during  the  imaging  time 
of  exposure,  degradation  results.  Examples  of  degra¬ 
dation-inducing  motions  are  random  and  determinis¬ 
tic  camera  vibrations,  camera  rotation,  and  linear 
camera  motion  along  a  flight  path. 

2.  Turbulent  media.  Differences  in  air  temperature  at 
varying  altitudes  cause  variations  in  the  air’s  index 
of  refraction  and  give  rise  to  the  random  phenome¬ 
non  turbulence.  Hufnagel  and  Stanley1  have  char¬ 
acterized  the  effect  of  this  type  of  degradation  for 
both  long-  and  short-term  exposures. 

3.  Diffraction-limited  optics.  The  optical  elements  in 
the  imaging  system  have  a  filtering  effect  on  the 
resulting  image.2  Under  very  low  noise  condi¬ 
tions,  even  the  filtering  effect  caused  by  an  ideal 
imaging  system  may  be  reduced  by  restoration  and 
resolution  beyond  the  diffraction  limit  ob¬ 
tained.'  5 

4.  Optical  aberrations.  Aberrations  present  in  the  im¬ 
aging  system  also  degrade  the  output  image.  Resto- 


‘R.  E.  Hufnagel  and  N.  R.  Stanley,  “Modulation  Transfer 
Function  Associated  with  Image  Transmission  through  Tur¬ 
bulent  Media,”  J.  Opt.  Soc.  Am.  54,  52-61  (1964). 

J.  W.  Goodman,  Introduction  to  Fourier  Optics,  McGraw- 
Hill,  San  Francisco  (1968). 

T.  S.  Huang,  W.  F.  Schreiber,  and  O.  J.  Tretiak,  “Image 
Processing,”  Proc.  IEEE  59,  1586-1609  (1971). 

JB.  R.  Frieden,  “Restoring  with  Maximum  Likelihood  and 
Maximum  Entropy,”  J.  Opt.  Soc.  Am.  62,  51 1-518  (1972). 

'B.  R.  Frieden  and  J.  .1.  Burke,  “Restoring  with  Maximum 
Entropy  II:  Superresolution  of  Photographs  of  Diffract'on- 
Blurred  Impulses,"  J.  Opt.  Soc.  Am.  62,  1202-1210(1972). 
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ration  techniques  have  been  developed  that  reduce 
the  effects  of  defocus,  coma,  curvature-of-field, 
astigmatism,  and  distortion  aberrations.1'' 

5.  Thermal  effects.  Temperature  differentials  across 
the  imaging  system  elements  and  overall  operating 
temperature  changes  can  cause  image  degradation. 
Spherical  aberration  and  a  shift  in  focus  have  been 
shown  to  arise  from  thermal  effects  caused  by 
aerodynamic  heating  in  infrared  missile  seeker  sys¬ 
tems.'' 

These  are  only  a  feu  of  the  more  important  degrad¬ 
ing  effects.  Additional  sources  of  degradation  include 
chromatic  aberration,  detector  nonuniformities,  camera 
shutter  effects,  and  scan  jitter  in  video  systems. 

1.3  APPLICATIONS  OF  RESTORATION 

Image  restoration  can  be  applied  to  many  interest¬ 
ing  image  problems  in  a  wide  variety  of  fields.  Most 
applications  today  are  restricted  to  data  that  have  been 
stored  on  magnetic  tape  or  some  other  storage  medium 


and  processed  long  after  the  image  was  formed.  How¬ 
ever,  the  development  of  rapid  algorithms  and  advanced 
hardware  implementations  should  make  real-time  image 
restoration  feasible. 

One  of  the  first  applications  of  image  restoration  was  j 

in  the  1960s  at  the  California  Institute  of  Technology 
Jet  Propulsion  Laboratory,  images  returned  from  the 
Mariner  spacecraft  contained  geometric  distortion  creat-  | 

ed  by  the  vidicon  on-board  camera.10  Digital  restora¬ 
tion  techniques  were  used  to  remove  the  distortion.  The 
algorithm  worked  by  locating  registration  reseau  marks  i 

and  calculating  a  coordinate  transformation  that  was  j 

subsequently  applied  to  the  image. 

Since  then,  restoration  techniques  have  been  applied 
in  the  diverse  areas  of  medicine  (X  rays,  acoustic  im¬ 
agery),  surveillance  data  (satellite  and  aircraft  imagery), 
oil  exploration  (seismic  signals),  and  forensic  science 
(smudged  fingerprints).  Restoration  has  even  found  ap¬ 
plication  in  the  music  world  as  Stockham  demonstrated  J 

by  restoring  Enrico  Caruso  recordings  using  homo¬ 
morphic  deconvolution." 

i 

1 

( 


2.0  IMAGING  SYSTEM  MODELS 


I  o  evaluate  and  reduce  the  effect'  ol  the  image  degra¬ 
dations  discussed  m  Section  1,0.  a  mathematical  model 
tot  the  imagine  stem  :s  needed.  Once  a  system  model 


(i  M  Kohl'iits  and  I  S  |  In, me.  "Imei'c  tillering  lot 
linear  Shut  \  atiant  hnaging  Svsicni'."  Pm,  HII  ML 
.'tv;  s-:  i  pc; i 

\  V  ''.re  Jink.  '  Spme  \  it  i.iiit  Imaji  Restoration  i'> 
<  oordmaic  I  rati'torm.iuon-."  /  ()[>!.  So<  \m  64.  I 's  14-1 
l  I T4  i 

\  \  s.iaJiu-  inJ  M  I  IVvinvi.m.  "Kesiot .limn  ol 

Vitetii, til'll!  uni  <  .it .it lit o  o!  I  lelsi."  /  (l/'t.  Soi  Ini.  65. 
"i;  'T  I  IT' i 

t  I  ilitiris.  ”1  valuation  , >1  f  nvironinentalOptie.il  f  Heel' 
>>!t  theli  Speevl  Heiin-rhet isal  Domes  I  sing  a  k,t\  I  race 
Model."  HU  M’l  llli:)»6.|  i if, 6  I  \pt  I9S6) 


lias  been  specified,  analytic  methods  can  be  used  to  de¬ 
rive  restoration  techniques.  In  this  section,  we  first  de¬ 
fine  the  imaging  system  in  general  terms  and  then  impose 
vat tous  lestrieuoiis  on  itte  general  system  description  that 
lead  to  analytically  tractable  models. 


I)  V  O'Mandlcv  and  \\ .  It.  Cireen.  ‘‘Recent  Developments 
in  Digital  linage  Processing  at  the  Image  Processing  Labo¬ 
ratory  at  the  let  Propulsion  I  aboratorv."  Proc.  ILIE  60, 
s;i  s;s  i it;». 

t  (,  Stockham.  T.  M.  Cannon,  and  R.  ft.  Ingebretscn, 
"Itlind  Deconvolution  through  Digital  Signal  Processing," 
Proc  11  1  (  63.  T*  64;  (1975). 
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2.1  CANONIC  IMAGING  SYSTEM  MODEL 

A  canonic  model  of  image  formation  and  detection  is 
shown  in  Fig,  I.12  The  detected  image  plus  noise  is 

g(x,,yl )  =  5,  lb(x,,y,)\  +  pi,  (x,,y,) 

+  S2[b(,x,,y, )  |  n2(xity,)  .  (U 

where  rt,  and  n2  represent  signal-independent  noise, 
/t,  =  S;  \h\n2  is  signal-dependent  noise,  S,  (  • )  is  the 
detector  response  function,  S2  { •  |  is  the  signal-depen- 
dent  noise  function,  and  denote  coordinates  in 
the  image  plane.  For  this  canonic  model,  the  functions 
S,  and  S2  may  be  nonlinear. 

The  image  radiant  energy,  b(x,,y, ),  is  found  by  op¬ 
erating  on  the  object  radiant  energy  by  the  system  point 
spread  function  (psf),  h,  which  will  be,  in  general,  a 
nonlinear  function  of  spatial,  temporal,  and  frequency 
coordinates,  as  well  as  the  object  function,/,  i.e. , 

b(x, X, )  = 

\  \  |  \  h{x, ,>•„/, ,\,x„,y„,t0,\0>  f(x0.y0,t0,\0)) 

iiit 

x 


where  ( xQ  ,y0 )  are  the  coordinates  in  the  object  plane. 
Although  imaging  systems  are  nonlinear  (e.g.,  typical 
X-ray  photographs),  many  are  approximately  linear 
when  within  some  operating  constraints.  To  obtain  a 
more  analytically  tractable  model,  we  make  the  sim¬ 
plifying  assumption  of  linearity  and  neglect  the  chro¬ 
matic  effects. 

If  we  restrict  the  temporal  effects  to  relative  motion 
between  the  imaging  camera  and  the  object  during  the 
time  of  exposure,  T,  they  can  be  incorporated  into  the 
spatial  coordinate  dependency.  For  example,  for  ob¬ 
ject  motion  described  by  the  functions 

x,  =  m,  (x 0,y0,t0)  ,  (3) 


y,  =  m2(x0,y0,t0)  , 


inverse  functions  can  be  found  such  that 

to  =  ki  (x„,y„,x,)  =  k2(x0,y0,y,)  ,  (4) 

and  the  image  radiant  energy  is  described  only  in  terms 
of  the  object  and  image  spatial  coordinates. 


x  dx„  dy„  dtn  dk„  , 


(2) 


n 2<*/.  V,) 


"!<*/.  Vl'l 


Figure  1  Canonic  imaging  system  model. 


H  C  .  Andrews  and  B.  R.  Hunt,  Digital  Image  Restoration, 
Prentice-Hall,  Englewood  Cliffs,  N..T.  (1977). 
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2.2  I 'TEAR  SPATIALLY 
V  ‘  "lANT/IN  VARIANT  MODELS 

The  simplifications  discussed  in  Section  2.1  result  in 
the  expression  for  a  linear  spatially  variant  (LSV)  im- 
aeina  system, 


h(x..y, )  =  \  \  h(x,,y,,x„,yu  ) 


x  /(  v„,y„  )  dx„  cly, 


where  h  is  a  function  of  all  four  spatial  variables  and 
has  been  modified  to  incorporate  any  relative  motion 
between  object  and  camera.  A  common  example  of  an 
LSV  imaging  system  is  one  where  the  camera  moves 
with  constant  linear  velocity  normal  to  the  optical  axis 
when  imaging  a  three-dimensional  scene.  Objects  in  the 
foreground  w  ill  be  blurted  more  than  objects  in  the  dis¬ 
tance.  resulting  in  spatially  variant  degradation. 

Further  simplification  results  if  the  system  obeys  the 
property  of  superposition.  Then  the  system  is  linear  spa¬ 
tially  invariant  (LSI)  and 

hi x  ,y  )  =  \  \  h(x  -  ,v„,  .v,  -  y„ ) 

x 

x  /  ( .v,  ,y„  )  r/.v„  dy„  .  »6) 


which  is  the  familiar  two-dimensional  convolution  in¬ 
tegral.  This  representation  is  particularly  attractive  since 
it  lends  itself  to  F  ourier  analysis.  Consequently,  many 
restoration  techniques  have  been  developed  for  LSI  im- 
acint!  svstems. 


2.3  SPATIALLY  SEPARABLE  MODELS 

l  Hdcr  certain  circumstances,  the  imaging  system  psf 
may  be  spatially  separable  and 


h(x„y„x0,y„)  =  hx(x,,xu)  h2(y,,y0)  (LSV) 

=  /?,  (a-,  -  xa )  h2(y,  -  yu)  (LSI). 


A  separable  LSI  psf  might  arise,  for  example,  from  a 
unit-square-aperture  ideal  imaging  system, 


h(x, ,y, ,x„ ,v„  ) 


r  sinirU,  -  x„)1 
L  v(x,  -  xa )  J 

X  fiinir(^  -^'1  .  ,8, 

L  rr(  (y,  -  ya)  J 


Separability  offers  the  advantage  of  performing  two- 
dimensional  restoration  by  sequential  one-dimensional 
operations.  This  results  in  a  tremendous  simplification 
in  the  computer  implementation  of  the  restoration  al¬ 
gorithm. 

Other  restrictions  may  be  put  on  the  imaging  system 
such  as  nonnegativity  to  ensure  positive  images  or  loss¬ 
less  imaging  constraints.  To  coincide  with  digitized  com¬ 
puter  array  representations,  we  also  need  to  develop  a 
discrete  system  model.  A  concise  discrete  representation 
lexicographically  orders  the  N  by  N  two-dimensional 
object,  noise,  and  image  data  into  one-dimensional  vec¬ 
tors  each  of  length  N2.  The  system  psf  can  then  be 
modeled  as  an  N':  by  /V:  matrix  of  values  [H\  result¬ 
ing  in 

b  =  [//]  f.  (9) 

The  matrix  \H]  can  be  shown  to  be  block  Toeplitz 
for  LSI  imaging  systems. i:  A  matrix  is  Toeplitz  if  the 
entries  on  each  diagonal  have  the  same  value  (i.e., 
h„  -  hkl  for  /  -  j  =  k  -  /);  a  matrix  is  block 
Toeplitz  if  it  has  a  Toeplitz  partitioning  and  each  par¬ 
tition  submatrix  is  also  Toeplitz.  A  summary  of  the 
resulting  imaging  models  (noise-free  case)  for  both  con¬ 
tinuous  and  discrete  systems  is  given  in  Table  1. 


z.  /.  /,  s  s.  f  j- 
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Table  1 

Summary  of  imaging  system  models. 

Continuous  system  Discrete  system 


Noise-free  imaging  g(x,,y,)  =  \\h(x„y„r0,y0  )f(x0,y0  )dx0dyo  gv  =  £  £  hi)klfk,  or 
model  *  ' 

g  =  [H]i  (lexicographically  ordered) 


Lossless  incoherent  g(x,,y,)  >  0 
imaging  f(x0,y„)  >  0 

h(x„y„x„,y0)  >  0 
Jj h(x„y„x0,y0)  dx,  dy,  =  1 


Linear  spatially 
variant 


Nonseparable 

h(Xj,y„x„,y0)  -  h(x„y„x0,yo) 
Separable 

h(x.,y,,x0,y,J )  =  h{  (x,,x„  )h2  (y,  ,y„  ) 


g,j  -  0 

fki  2:  0 

hij.ki  >  0 

X)  ]C  ^Ij.kl  =  1 

‘  J 

im  =  [hi 

[Hi  =  [//,]  ®[//2] 


Linear  spatially 
invariant 


Nonseparable 

h(x,,y,,x0,y„)  =  h(x,  -xu,y,  -y„  ) 
Separable 

h(x„y„x0,y0  )  =  h,  (x,  -x0  )h2  (y,  -y„  ) 


[H]  =  block  Toeplitz 

[H]  =  [//,]®[//2] 
[HA,  \H2  \  Toeplitz 


3.0  A  POSTERIORI  DEGRADATION  DETERMINATION 


Restoration  as  defined  in  Section  1 .0  requires  a  priori 
knowledge  of  the  degradation  performed  on  the  original 
object.  Most  of  the  restoration  techniques  described  in 
Section  4.0  require  knowledge  of  the  imaging  system 
psf.  In  practical  restoration  applications,  we  are  com¬ 
monly  faced  with  the  task  of  a  posteriori  determination 
of  the  degrading  psf  from  the  output  blurred  image. 


Various  approaches  for  accomplishing  restoration  are 
possible.  The  techniques  outlined  below  are  restricted 
to  LSI  imaging  system  applications.  However,  if  the  sys¬ 
tem  is  LSV,  these  techniques  can  be  extended  by  divid¬ 
ing  the  image  into  isoplanatic  patches  to  determine  the 
local  psf  in  each  region  (see  Section  4.3.1  for  further 
discussion).  If  the  original  image  has  a  specific  struc- 
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ture,  direct  measurement  of  the  psf  can  be  performed. 
Point,  line,  and  edge  sources  imaged  by  the  system  can 
be  used  to  directly  measure  or  calculate  the  system 
psf. 11  Alternately,  psf  estimation  can  be  accomplished 
by  homomorphic  signal  processing  techniques  per¬ 
formed  on  the  blurred  image.  If  a  specific  parametric 
form  for  the  degradation  is  assumed  (e.g.,  linear  con¬ 
stant  velocity  blur,  simple  optical  defocus),  estimation 
techniques  that  determine  the  model  parameters  can  be 
used  to  develop  an  expression  for  the  system  psf. i: 

3.1  DEGRADATION  PSF  MEASUREMENT 
TECHNIQUES 

There  are  three  direct  measurement  techniques  that 
can  frequently  be  used  to  determine  the  system  psf. 

3.1.1  Point  Source  Measurement 

Locate  an  image  of  an  isolated  point  source  in  the 
degraded  output  and  measure  the  psf  directly  by 

x 

g(.v  ,.t-  )  =  \  \  h(x,  -  x„,  v,  -  v„ ) 

-  X 

X  f>(x„,y„)  dx„  dy„  =  h(x,,\\)  .  (10) 

This  method  is  particularly  applicable  in  astronomical 
imaging  with  stellar  point  sources. 

3.1.2  Fine  Source  Measurement 

Locate  an  image  of  an  isolated  line  source  in  the 
blurred  output  and  measure  the  line  spread  function  (lsf) 
directly  by 

x 

L'(.v  ,.r  )  -  \  \  h(x,  -  x„,  y,  -  v„ )  5(x,) 
dx„  dy„ 

■  x 

=  I  /t(.v,,v,  -  i1,,  I  dv„  =  h,( )  . 

•  -  '  '  '  (1!) 

In  the  Fourier  domain, 

/T|.g(.v. ,  v,  ) )  =  G(J\ )  =  HU\. 0)  .  (12) 

where  FT  1  •  |  denotes  a  two-dimensional  Fourier  trans¬ 
form.  If  the  psf  is  known  to  be  circularly  symmetric 

V  Rosen feld  and  A.  C.  Kak,  Digital  Picture  Processing. 
2nd  ed.,  -Xeadcniie  Press,  New  York  (19X2). 


(an  isotropic  function),  H  and  therefore  h  are  deter¬ 
mined.  Otherwise  we  must  Find  the  lsf  for  line  sources 
positioned  at  varying  angles  in  the  image  to  determine 
H  along  several  radial  lines.  Interpolation  is  used  to  ob¬ 
tain  the  transfer  function  values  at  rectangular  grid  lo¬ 
cations. 

3.1.3  Edge  Source  Measurement 

Locate  an  isolated  edge  in  the  degraded  image  and 
measure  the  edge  spread  function  (esf).  The  derivative 
of  the  esf  is  the  lsf  and,  therefore, 

H/  (/,.,/.. ) 

He  WxJv)  =  —  7  (edge  along  y  axis)  .  (13) 

i2  *fx 

From  this  result,  we  can  determine  the  psf  as  before. 
Par-target  test  patterns  imaged  by  the  system  are  fre¬ 
quently  used  for  this  technique. 


3.2  DEGRADATION  PSF  ESTIMATION 
TECHNIQUES 

The  techniques  described  below  are  based  on  homo¬ 
morphic  signal  processing  (HSP)  and  cepstral-like  tech¬ 
niques  to  estimate  the  degrading  filter  and  its  param¬ 
eters. 

3.2.1  HSP  Filter  Estimation 

The  blurred  image  is  subdivided  into  N  regions.  This 
technique  requires  that  the  nonzero  extent  of  the  psf 
be  small  compared  to  the  area  of  each  region.  Each  re¬ 
gion  is  Fourier  transformed  to  yield 

G,  (/,,/,)  a  HU\.L)  F,  (fxJy  ).  (14) 

where  i  =  1,  2,  ....  N.  The  logarithm  of  the  magni¬ 
tude  is  taken  and  the  result  is  summed  over  the  N  regions 
by 

V  \ 

2>  |G,|  =  N  In  \H\  +  V/n  \F,\  .  (15) 

I-  1  /  -  1 

The  degrading  filter  is  estimated  as 

I"!  =  exp[^  ln\G,\  -  ][>  |F,  |  ] 


f  ./ 

s' 
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For  this  technique  to  work,  the  object  Fourier  trans¬ 
form  in  each  region  must  either  be  known  or  estimated. 
Another  possibility  is  to  assume  that  for  large  enough  A/, 


by  the  linear  constant  velocity  motion  of  the  camera 
during  exposure.  Equations  3  and  4  (or  others  of  that 
form)  can  be  used  to  determine  the  LSI  psf, 


'£tln\Fl  |  s  Constant. 


The  conditions  that  must  be  met  to  satisfy  this  last 
criterion  are  not  well  defined  and  require  further  inves¬ 
tigation. 

3.2.2  Power  Spectral  Density  Estimation 

A  related  approach  treats  the  object  and  image  as 
stochastic  processes  and  uses  the  function  power  spec¬ 
tral  densities  (psd)  to  determine  the  degrading  Filter.  For 
linear  systems, 

S„  =  \H\2  S„  ,  (18) 

where  Seit  and  Sff  are  the  psd  of  the  image  and  object 
functions,  respectively.  The  blurred  image  psd  can  be 
estimated  by  the  spatial  average 

(19) 

■'\=l 

where  the  assumption  is  made  that  the  stochastic  im¬ 
age  process  is  ergodic.  The  resulting  filter  is 


;,Dc'iT  = 


The  object  psd,  Srl ,  is  either  assumed  known  or  is  also 
estimated  as  a  spatial  average,  and  the  IF,!, 
i  —  1,  2,  .  .  ,  A  are  assumed  known. 

Both  of  the  above  estimation  techniques  lead  to  an 
expression  for  the  magnitude  of  the  degrading  filter; 
they  give  no  information  about  the  filter  phase.  For 
degrading  filters  that  are  known  a  priori  to  be  real  and 
to  contain  no  zero  crossings  (e.g.,  Gaussian  blur),  the 
filter  has  been  sufficiently  determined.  If  the  filter  trans¬ 
fer  function  does  contain  zero  crossings,  we  may  be  able 
to  postulate  a  parametric  form  for  the  filter.  Fstimation 
techniques  can  then  be  used  to  determine  the  parame¬ 
ters.  Two  examples  of  this  parametric  estimation  meth¬ 
od  arc  given  below. 

3.2.3  Parametric  Filter  Estimation 

3.2.3. 1  Linear  constant  velocity  blur.  We  consider 
an  imaging  system  where  the  only  degradation  is  caused 


h(x,y)  = 


-  0  <  xcos<t>  +  ysin<£  <  vT 
v  and  y  =  xtan<t> 

0  elsewhere. 


The  filter  transfer  function  is 


H{fxJy)  =  T 


sin(ir  fvT)e  nf'T 
Tr/v'r 


where  T  is  the  time  of  exposure,  v  is  the  linear  velocity, 
<t>  is  the  angle  of  motion  relative  to  the  x  axis,  and 
/  =  fx  cos  <t>  +  fy  sin  0. 

This  filter  clearly  has  zero  crossings  and  a  phase  func¬ 
tion  associated  with  it.  The  HSP  techniques  discussed 
in  Section  3.2  will  not  work  for  this  degradation  and 
some  other  method  of  estimating  the  filter  is  required. 
In  particular,  we  need  to  estimate  the  parameters  vT 
and  4 >.  To  accomplish  this,  we  compute  one-dimensional 
Fourier  transforms  of  the  log  magnitude  of  the  blurred 
image  spectrum  with  respect  to  both  the  fx  and  fy  fre¬ 
quency  axes.  The  periodic  structure  of  the  filter  transfer 
function  gives  rise  to  peaks  in  the  two  one-dimensional 
Fourier  transforms  thus  generated  at  the  values 
x  -  vT  sin  4>  (Fourier  transform  with  respect  to 
and  y  =  vT  cos  4>  (Fourier  transform  with  respect  ti 
/, ).  This  cepstral-like  technique  allows  us  to  determine 
both  unknown  parameters,  vT  and  4>,  which  in  turn  uni¬ 
quely  specify  the  degrading  filter. 

3.2.3. 2  Simple  defocus  aberration.  The  degrada¬ 
tion  arising  from  a  badly  defocussed  optical  system  can 
be  modeled  as  the  Airy  disk  pattern12 


H{fX.fy)  = 


A  4}  +7?' 


where  7,  is  a  Bessel  function  of  the  first  kind  and  A 
is  the  effective  radius  of  the  circular  aperture.  The  log 
magnitude  of  the  image  spectrum  after  inversion  and 
clipping  yields  large  positive  values  forming  a  circle  of 
radius  A.  With  this  technique,  the  radius  can  thus  be 
estimated  and  the  filter  determined.  Obviously,  care 
must  be  taken  with  both  parametric  estimation  methods 
to  avoid  logarithms  of  zero  values. 
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3.3  NOISE  CONSIDERATIONS 

So  far,  no  mention  of  the  effect  of  noise  has  been 
made,  in  general,  the  model  for  the  output  image  will 
consist  of  the  linearly  filtered  object  plus  signal  inde¬ 
pendent  noise, 

K  -  h  *  f  +  n 

and 

G  =  HF  +  V  .  (24) 

The  introduction  of  noise  into  the  system  will  clearly 
degrade  the  measured  and  estimated  values  for  the  sys¬ 
tem  filter  just  derived. 

One  way  to  reduce  the  effect  of  noise  is  multiframe 
averaging.  Several  frames  of  the  degraded  image  may 
be  available  (as  is  often  the  case  with  video  data),  and 
noise  reduction  can  occur  by  averaging  these  image 
frames.  Multi  frame  averaging  relics  on  the  assumption 


that  the  detector  noise  is  independent  from  frame  to 
frame  to  effect  a  reduction  in  the  noise  variance. 

The  imaging  system  psf  can  also  be  considered  to  be 
varying  from  frame  to  frame  while  the  object  remains 
unchanged.  This  is  the  case  when  detector  camera  mo¬ 
tion  or  random  turbulent  media  variations  are  incor¬ 
porated  in  the  psf.  The  summation  over  J  image  frames 
yields 

g'  =  h'  *  f  +  n'  ,  (25) 

where  the  new  noise  term,  n' ,  retains  the  same  mean 
as  the  original  noise  but  has  a  variance  that  is  reduced 
by  the  factor  \/J.  The  new  system  psf  is  h'  and  is  the 
ensemble  average  of  the  individual  frame  psf.  If  the 
noise  n'  has  a  known  mean  value  and  the  variance  is 
sufficiently  small,  degradation  measurement  or  estima¬ 
tion  techniques  can  be  applied  to  the  averaged  image 
in  Eq.  25  to  generate  an  estimate  of  h' .  Restoration 
can  then  be  carried  out  on  this  averaged  image  to  ob¬ 
tain  a  restored  version  of  the  original  object. 


4.0  IMAGE  RESTORATION  TECHNIQUES 


When  surveying  the  image  restoration  literature,  it 
becomes  apparent  that  a  plethora  of  restoration  tech¬ 
niques  exists.  Each  technique  is  based  on  particular  as¬ 
sumptions  about  the  object  and  imaging  system  and  is 
motiv  ated  by  some  form  of  restoration  criteria.  A  com¬ 
prehensive  discussion  of  all  restoration  methods  is  be¬ 
yond  the  scope  of  this  report.  Instead,  we  consider  only 
the  most  popular  and  successful  methods  for  both  ESI 
and  I.SV  imaging  systems.  The  strengths  and  weaknesses 
of  each  method  are  identified  and  restoration  examples 
are  supplied  where  possible;  however,  rather  than  pre¬ 
senting  a  somewhat  long  and  disjointed  list  of  restora¬ 
tion  techniques,  we  begin  with  the  general  concepts  that 
tie  these  varied  methods  together. 

4.1  RESTORATION  CRITERIA 

(liven  a  degraded  image,  we  would  like  to  use  the 
prior  know  ledge  obtained  about  the  imaging  system  to 


perform  some  operation  on  the  degraded  image  to  ef¬ 
fect  an  “improvement”;  i.e.,  a  restored  image  is  gener¬ 
ated  that  in  some  sense  more  closely  resembles  the  origi¬ 
nal  object.  The  choices  of  criteria  used  to  motivate  the 
improvement  and  determine  the  restoration  technique 
fall  into  four  general  categories:  (a)  least  squares  esti¬ 
mates,  (b)  equivalent  power  spectral  densities,  (c)  Bayes¬ 
ian  estimates,  and  (d)  ad  hoc  methods. 

Techniques  based  on  least  squares  estimation  treat 
the  original  object  as  a  deterministic  unknown  function. 
This  criterion  requires  that  the  squared  error  between 
the  degraded  image  and  the  filter  d  estimate  be 
minimised, 

Min|e j  =  Min  \  \g  -  h  *  ,f\:\  .  (26) 

/'  i 

Because  of  the  invariant  filtering  operation,  least  squares 
estimation  only  results  in  ESI  restoration  techniques. 


16 


AjV.jV'j'W.v.V-  •.  ,N  , 


tfiwiliaift 


THE  JOHNS  HOPKINS  UNIVERSITY 

APPLIED  PHYSICS  LABORATORY 

LAUREL.  MARYLAND 


The  second  category  is  that  of  equivalent  power  spec¬ 
tral  densities.  Techniques  that  fall  into  this  category  treat 
the  object  as  a  stochastic  process  and  generate  a  restored 
image  whose  psd  is  equal  to  the  known  original  object 
psd, 

%  -  Sff  •  (27) 

Bayesian  estimation  gives  rise  to  a  number  of  estima¬ 
tion  criteria  depending  on  the  cost  function  assigned  to 
errors. 14  The  stochastic  properties  of  the  object  and 
noise  are  used  to  develop  minimum  mean  squared  error 
(MMSE),  maximum  entropy,  and  maximum  a  posteriori 
(MAP)  estimates. 

In  addition  to  techniques  based  on  these  well-defined 
criteria,  others  developed  on  heuristic  arguments  can 
also  provide  good  restoration.  Application-specific  con¬ 
straints  such  as  non-negativity  or  second  order  smooth¬ 
ness  can  also  be  introduced  to  generate  various  hybrid 


techniques.  Table  2  lists  the  categories  of  restoration 
criteria  and  the  restoration  techniques  associated  with 
each  one.  These  techniques  are  discussed  in  detail  in 
Sections  4.2  and  4.3. 


4.2  LINEAR  SPATIALLY  INVARIANT 
TECHNIQUES 

Recall  that  the  LSI  imaging  system  model  is 


g(x,,y,)  =  j  J  h(x,  -  x0,  y ,  -  ya) 

—  00 

x  f(x0,  y0  )  dx0  dy0 
+  n{x,,y,)  ,  (28) 


Table  2 

Restoration  criteria  and  techniques. 


Category 


Criterion 


Techniques 

LSI  LSV 


1.  Least  squares  estimate 


Mini  |g  -  h  *  /| 2 1  Inverse  filtering 


2.  Equivalent  power  spectral 
density 


%  ~  Sfi 


Homomorphic  filtering 
Geometric  mean  filtering 
(a  =  1/2,  7  =  1) 


3.  Bayesian  estimate 

Minimum  mean  square  Min  E\  \f  -  f\ 1 1 
error 


Wiener  filtering  Analytic  continuation 

Recursive  filtering 


Maximum  entropy 


Max  |  -  /  In  f\ 


Maximum  entropy 
filtering 


Maximum  a  posteriori  Max  Pr\f\g\ 


Maximum  a  posteriori 
filtering 


4.  Ad  hoc  methods  Constrained 

deconvolution 
Geometric  mean  filtering 

(0  s  a  <  1,7) 


H.  I  Van  Trees,  Detection,  Estimation,  and  Modulation 
Theory,  Part  I,  Wiley  &  Sons.  New  ''  ’  (1968). 
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r 


where  the  noise  is  modeled  as  additive  and  signal  inde¬ 
pendent.  In  the  Fourier  domain,  this  becomes 


GCA./,)  =  FV\J„  )  H(fxJy)  +  N(fxJy)  .  (29) 


The  subsections  below  discuss  various  methods  of  im¬ 
age  restoration  applicable  to  LSI  systems.  They  incor¬ 
porate  techniques  that  treat  the  degrading  filter,  H ,  as 
deterministic  and  known,  deterministic  with  unknown 
parameters,  and  completely  random. 

4.2.1  Inverse  Filtering 

The  technique  of  inverse  filtering  or  method  of  least 
squares  assumes  a  known  (or  accurately  estimated)  de¬ 
terministic  filter  transfer  function,  H.  It  generates  a  lin¬ 
ear  restoration  filter  with  psf  m(x,y),  which  satisfies 
the  criterion  of  minimum  squared  error, 

Min| el  =  Min  I  |g  -  h  *  f\ 2 )  .  (30) 


f  -  m  *  g  . 


This  is  easily  solved  to  yield* 


V/lllu.rsc  =  FT\  m  I  =  \/H 


/W  =  FT\f\  =  G/H  .  (33) 

Equation  33  demonstrates  that  the  inverse  filter  has 
severe  problems  whenever  H  =  0.  If  H  has  zeros  in  the 
desired  image  bandwidth,  the  image  cannot  be  perfectly 
restored,  even  in  the  absence  of  noise,  because  of  the 
indeterminate  0/0  ratios  that  occur.  When  noise  is  pres¬ 
ent,  the  zeros  of  H  serve  to  amplify  the  noise  power 
at  these  frequencies.  If  H  is  reasonably  band  limited, 
high-frequency  noise  power  can  also  be  severely  in¬ 
creased 

A  rough  measure  of  the  degradation  in  signal-to-noise 
ratio  (SNR)  from  pre-  to  postrestoration  is  found  by 
defining  a  voltage  SNR  for  the  degraded  image, 


S\R,  = 


!  \h  *  /|| 


where  the  energy  conserving  property  of  h  has  been 
used12  and  1 1  •  1 1  denotes  the  signal  norm.15  The  SNR 
for  the  image  restored  by  inverse  filtering  is  defined  as 


SNRf  = 


| m  *  n\ 


Nil  IN 


TT  •  (35) 


The  ratio  of  the  two  SNRs  is  the  degradation  in  SNR 
caused  by  the  inverse  filtering  restoration  process, 

SNRf  1  1 

2  =  77  TT  ®  7^rr^7  ,  (36) 


where  Parseval’s  theorem  accounts  for  the  last  equali¬ 
ty.  Typical  values  for  ||//_1||  for  real  imaging  sys¬ 
tems  are  on  the  order  of  100  or  higher,  demonstrating 
that  the  SNR  can  be  severely  degraded  by  a  factor  of 
100  or  more  when  using  inverse  filtering  restoration. 

This  points  out  the  general  ill-conditioned  nature  of 
image  restoration,  which  results  from  the  property  that 
small  perturbations  in  the  degraded  image,  g(x,y),  can 
cause  large  changes  in  the  restored  image,  f(x,y).  For 
example,  in  the  noise-free  case, 

g  =  ft  *  f 

f  =  m  *  g  =  /  (assuming  no  0/0  ratios),  (37) 
but  when  an  arbitrarily  small  amount  of  noise  is  added, 
g  =  h  *  f  +  n, 

f=m*g=fy-nh  ,  (38) 

where  n6  is  not  necessarily  small  and  can  in  fact  be 
quite  large. 

Inverse  filtering,  however,  does  have  its  advantages. 
It  is  easy  to  implement  and  can  be  done  quickly.  It  re¬ 
quires  only  knowledge  of  the  system  psf,  unlike  many 
other  techniques  that  require  knowledge  of  the  noise 
and  object  characteristics.  In  addition,  in  high-SNR  en¬ 
vironments  it  gives  restorations  with  good  resolution 
(providing  there  are  few  zeros  of  H  in  the  image  band¬ 
width).  To  avoid  indeterminate  ratios,  a  pseudoinverse 
filter  caii  be  defined  as12 


1  Pscudotn  verse 


=  H  =  lim 


y-o  \H\l  +  y 


>‘L.  E.  Franks,  Signal  Theory,  rev.  ed.,  Dowden  &  Culver. 
Stroudsburg,  Pa.  (1981). 
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Figure  2  shows  an  original  object  function.  Figure  3 
shows  the  object  after  it  has  been  blurred  by  a  low-pass 
filter  and  degraded  by  noise  for  both  a  high  SNR 
(33  dB)  and  a  lower  SNR  (23  dB).  The  restoration  gener¬ 
ated  by  inverse  filtering  is  given  in  Fig.  4  for  the  high- 
and  low-SNR  cases.  As  expected,  the  amplified  noise 
severely  degrades  the  restoration  at  low  image  SNRs. 

4.2.2  Wiener  Filtering 

The  poor  noise  performance  of  inverse  filtering  led 
to  the  development  of  alternate  restoration  techniques 
designed  to  restore  degraded  images  with  lower  SNRs. 
This  approach  treats  the  original  object  and  noise  as 
statistically  uncorrelated  random  functions  and  con¬ 
structs  a  Bayesian  estimate  of  the  object.  The  system 
transfer  function,  //,  is  initially  treated  as  determinis¬ 
tic  and  known.  This  method  entails  finding  the  LSI 
restoring  filter  that  minimizes  the  mean  squared  error 
of  the  resulting  estimate.  The  linear  filter  that  accom¬ 
plishes  this  is  commonly  known  as  the  Wiener  filter. 16 

The  W'iener  filter  is  derived  as  follows.  The  mean 
squared  estimate  error  is 


e  =  E\  |/  -  /|  - 1 

=  f|  1/  -  rn  *  J?r  I  . 


’ '  %  v  V  ■ 


Figure  3  Object  with  Gaussian  blur  and  additive  noise 
with  (a)  SNR  =  33  dB  and  (b)  SNR  =  23  dB.12 

where  E\  ■ )  denotes  the  ensemble  average  and  m  is  the 
psf  of  the  Wiener  restoration  filter.  The  technique  re¬ 
quires  that  we  find  the  function  m(x,y),  which 
minimizes  this  error.  Using  the  central  concept  of  lin¬ 
ear  mean  square  estimation  theory,  i.e.,  the  orthogonal¬ 
ity  principle,  which  states  that  the  error  in  the  estimate 
must  be  orthogonal  to  the  data,  wc  obtain 


Figure  2  Original  object  function.'2 


E\  [fix,  ,.V|  )  -  m(.v,,V|  )  *  £(.V|,y,  )] 

£*(*,,>•;)  |  =  0  (41) 


A.  Papoulis,  Probability,  Random  Variables,  and  Stochas¬ 
tic  Processes,  McGraw-Hill,  New  York  (l%5). 


R,  (Jlv.Av)  =  m(\x,\v)  *  R  (Ay, Ay) 
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(a) 


(b) 


Figure  4  Inverse  filtering  restoration  of  Fig.  3  with  (a) 
SNR  =  33  dB  and  (b)  SNR  =  23  dB.’2 

where  R,„  and  RK„  are  the  cross-  and  autocorrelation 
functions  of/ with  g  and  g,  respectively.  The  Fourier 
transform  of  F.q.  42  yields 

■S',,  =  M  S„  .  (43) 


^Wiener  -  $fg  I $gg  •  (44) 

We  have  implicitly  made  the  assumption  that  the  ran¬ 
dom  object  and  noise  processes  are  statistically  wide- 
sense  stationary.  This  is  frequently  not  true  for  real  ob¬ 
jects,  thereby  limiting  the  class  of  objects  for  which 
Wiener  filtering  truly  results  in  the  MMSE  estimate. 
However,  we  proceed  by  not  only  assuming  that  f(x,y) 
and  n(x,y)  are  wide-sense  stationary,  but  also  constrain 
them  to  be  uncorrelated.  Then 


Mvicner 


\H\ 2  Sff  +  Sn„ 


(45) 


where  Snn  is  the  psd  of  the  noise.  This  may  also  be 
written  as  the  product  of  a  simple  inverse  filter  and  a 
modifying  filter, 


icner 


Mo 

H 


1 

H 


[r 


+  S„„/|//|2  S, 


'//- 


•  (46) 


The  modifying  filter,  Mot  is  dependent  on  the 
stochastic  properties  of  the  random  object  and  noise 
functions.  At  frequencies  where  the  noise  psd  is  negligi¬ 
ble,  the  Wiener  filter  behaves  as  an  inverse  filter.  At 
frequencies  where  the  noise  dominates,  the  modifying 
filter  provides  a  weighting  factor  that  adjusts  the  value 
of  M  to  be  appropriately  small.  The  modifying  filter 
thus  supplies  a  smooth  transition  between  the  two  noise 
extremes.  Slepian17  considers  the  case  where  the  de¬ 
grading  filter,  H,  is  also  modeled  as  a  random  function, 
and  finds  that  the  restoring  filter  is  found  by  replacing 
H *  and  \H\l  in  Eq.  45  by  £(//*)  and  E\\H\2\, 
respectively. 

The  Wiener  filter  does  not  exhibit  the  singularity 
problems  associated  with  the  inverse  filter  at  the  zeros 
of  the  system  filter  H.  In  fact,  it  is  the  presence  of  noise 
that  ensures  that  the  Wiener  restoration  filter  is  zero 
whenever  H  is  zero. 

The  restoring  filter  derived  in  Eq.  46  is  the  optimum 
LSI  MMSE  filter.  However,  nonlinear  filters  may  exist 
that  yield  smaller  mean  squared  errors.  In  general,  the 
MMSE  estimate  is  given  by  /  =  Ef/lg),  which  may 
be  a  nonlinear  function  of  g.  If  /  and  n  are  jointly 
Gaussian  and  stationary,  the  MMSE  estimate  reduces 
to  a  simple  linear  filtering  operation  and  the  true  MMSE 
estimate  is  equal  to  the  LSI  MMSE  estimate. 


where  S,..  and  .S’.,,,  arc  the  analogous  psd.  The  MMSE  1  D.  Slepian,  “Linear  Least-Squares  Filtering  of  Distorted  Im- 
Wiener  restoration  filter  is  ages."  ./.  Opt.  Soc.  Am.  57,  918-922  (1967). 
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Examples  of  Wiener  filtering  restoration  of  the 
degraded  images  in  Fig.  3  are  shown  in  Fig.  5.  Some 
resolution  is  lost  in  the  high-SNR  case;  however,  sig¬ 
nificant  improvement  is  made  in  the  low-SNR  case  (cf. 
Fig.  4). 

A  paramenic  version  of  the  MMSE  Wiener  filter  has 
also  been  developed.  For  a  specific  application,  the  user 
can  weight  the  psd  ratio  in  the  filter  demoninator  (Eq. 


* 


m 

in  in 


tiyir. 


Figure  5  Wiener  filtering  restoration  of  Fig.  3  with  (a) 
SNR  =  33  dB  and  (b)  SNR  =  23  dB.1J 


46)  to  reflect  the  degree  of  importance  that  should  be 
attached  to  this  stochastic  information.  This  results  in 
the  parametric  Wiener  filter 


'Parametric  Wiener 


■Mr 


+  yS„„  / 1  H\  2S,, 


where  7  >  0.  The  parameter  7  is  determined  subjec¬ 
tively  by  the  user  for  a  particular  image. 

4.2.3  Geometric  Mean  Filtering 

We  have  determined  that  the  simple  inverse  filter  can 
demonstrate  good  resolution  performance  at  spatial  fre¬ 
quencies  where  the  signal  dominates  the  noise  (usually 
lower  frequencies),  but  is  notoriously  poor  in  high-noise 
regions  (usually  higher  frequencies).  The  Wiener  filter 
has  excellent  noise  performance  but  achieves  this  at  the 
expense  of  smoothing  the  restored  image  by  the  modify¬ 
ing  filter,  Ma.  This  can  be  seen  by  examining  the 
MMSE  restored  image  in  the  Fourier  domain, 


=  4/wi,re,  G 
=  Ma  ( HF  +  AO 

H  NMg 

=  M„F  + - -  . 

u  v  r 


The  restored  image  is  no  longer  a  noisy  version  of  the 
ideal  object  as  in  inverse  filtering  restoration,  but  a 
smoothed  noisy  version  of  f(x,y ). 

Many  applications  require  better  combined  noise  and 
resolution  performance  than  either  method  provides  in¬ 
dividually.  A  heuristic  technique  that  attempts  to  recon¬ 
cile  the  trade-off  between  resolution  and  noise  is  geo¬ 
metric  mean  filtering  restoration.  The  restoration  fil¬ 
ter  is  defined  as 


^Geomcnic  mean  {^Inverse  1  [^Parametric  Wiener  1 

(49) 


The  parameters  0  <  a  <  1  and  7  >  0  are  user  defined 
for  a  particular  image.  The  user’s  choice  of  a  thus  shifts 
the  emphasis  from  Wiener  to  inverse  filtering  as  a  varies 
from  zero  to  one. 

4.2.4  Constrained  Deconvolution 

Both  Wiener  and  geometric  mean  filtering  require 
knowledge  of  the  stochastic  properties  of  the  noise  and 
the  original  object.  Constrained  deconvolution  is  simi- 


r*tr*,* 
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lar  to  inverse  filtering  in  that  it  treats  the  object  as  an 
unknown  deterministic  function. 

This  technique  minimizes  the  square  of  a  linear  con¬ 
straint  on  the  estimate, 

Min  f| c(x,y)  *  f(x,y)  |2|  ,  (50) 

subject  to  a  fixed  value  for  the  estimate  error, 

e  =  \g  -  h'f)2  =  E[n2\  .  (51) 


The  sampled  constraint  function,  £■(*,.>'),  is  frequently 
chosen  as  a  second  order  or  higher  difference  matrix. 
The  second  order  difference  matrix,  for  example,  is  de¬ 
fined  as  the  Kronecker  product  of  a  tridiagonal  matrix 
C,  with  itself,  C  =  C,  x  Cj .  The  matrix  C,  per¬ 
forms  the  one-dimensional  second  order  differencing 
operation  [/(.v+1)  -  /(*)]  -  l/(x)  -  f(x- 1)]  on 
the  data  and  is 


C, 


-2  1 
1  -2  1 

1  -2 


0 


0 


1 

-2  1 
1  -2 


Minimizing  the  second  order  difference  will  prevent  the 
resulting  estimate  from  containing  wild  oscillations. u  11 

The  method  of  Lagrange  multipliers  may  be  applied 
to  yield  the  restoration  filter 


Mi 


onstrained  deconvolution 


\H\ 2  +  >|C| 2 

i  r _ 1 _ 

H  Ll  +  7|C|2/|//|2 


where  C  =  FT\c\  and  y  is  related  to  the  Lagrange  mul¬ 
tiplier.  The  value  of  y  is  adjusted  so  that  the  fixed  er¬ 
ror  criterion  is  satisfied. 

Using  constrained  deconvolution,  one  can  thus  choose 
an  estimate  error,  E\n2 1  (perhaps  determined  a  poste¬ 
riori  from  the  degraded  image),  and  generate  a  restored 
image  that  has  the  desired  properties  provided  by  the 
constraint,  c.  The  attractiveness  of  this  technique  is  that 


the  restored  image  will  have  these  properties  but  they 
will  be  obtained  without  detailed  knowledge  of  the  ob¬ 
ject  and  noise  functions’  stochastic  characteristics. 

Figure  6  shows  the  image  restoration  achieved  by  each 
of  the  LSI  techniques  discussed  so  far  for  a  defocused 
object  with  additive  high-frequency  noise. 

4.2.5  Homomorphic  Deconvolution 

Homomorphic  image  processing  for  use  in  a  posteri¬ 
ori  degradation  determination  was  introduced  in  Section 
3.2.  We  can  proceed  with  the  results  of  that  section  to 
develop  a  homomorphic  restoration  filter.  This  tech¬ 
nique  is  the  work  of  Stockham  et  al. 11  and  Cannon,18 


(•)  <») 


Figure  6  Restoration  by  various  techniques;  (a)  origi¬ 
nal  object,  (b)  original  with  defocus  and  high-frequency 
noise,  (c)  inverse  filtering  restoration,  (d)  Wiener  filter¬ 
ing  restoration,  (e)  geometric  mean  filtering  restoration 
„  =  i/2,  -y  =  1,  (f)  constrained  deconvolution  restora¬ 
tion  using  C  =  second  difference.12 


"*M.  Cannon,  “Blind  Deconvolution  of  Spatially  Invariant 
Blurs  with  Phase,”  IEEE  Trans.  Acoust.  Speech  Signal  Pro¬ 
cess.  24,  58-63  (1976). 
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and  is  also  known  as  blind  deconvolution.  It  is  unique 
among  the  restoration  techniques  discussed  so  far  be¬ 
cause  it  does  not  assume  that  the  system  degradation 
filter  H  is  known. 

As  in  Section  3.2,  we  subdivide  the  degraded  image 
into  N  regions  and  Fourier  transform  each  region.  The 
psd  of  the  blurred  image  is  estimated  by  (ergodic  as¬ 
sumption) 


s»  a  -N  Eic.i1 


The  image  estimate  in  the  Fourier  domain  is 


F  =  MG  , 


where  M  is  the  homomorphic  restoration  filter  yet  to 
be  determined.  The  psd  of  the  estimate  is  therefore 


%  =  m 2  s„ 


The  restoration  criterion  (see  Table  2)  for  this  tech¬ 
nique  is  to  equate  the  original  object  and  estimate  im¬ 
age  psd,  i.e., 

%  =  Sff 


\W%  =  S/f  , 


which  yields 


'Homomorphic  deconvolution  I 


[ V««]W  •  (57) 


Here,  Sff  is  assumed  known  or  may  be  estimated  from 
an  undegraded  image  with  structure  similar  to  that  of 
the  original  object. 

The  homomorphic  restoration  filter  derived  is  the 
reciprocal  of  the  degrading  filter  estimated  by  Eq.  20. 
However,  it  is  not  a  simple  inverse  filter  as  the  follow¬ 
ing  argument  illustrates.  If  the  object  and  noise  are  wide- 
sense  stationary  and  uncorrelated  random  processes,  we 
know 

S„  =  \H\2  Sff  +  Sm  ,  (58) 


where  the  noise  is  assumed  to  be  additive  with  psd  Sm . 
Then, 


’  1  Urn  2sf/  +  snj  ■ 

which  is  the  magnitude  of  the  geometric  mean  filter  with 
a  =  1/2  and  7=1, 


I  ^Homomorphic  deconvolution  I  — 


,  a  =  1/2.  y  =  1  I  •  (60) 


'Geometric  mean. 


Thus  the  magnitude  of  the  homomorphic  filter  is  the 
magnitude  of  the  geometric  mean  of  the  inverse  filter 
and  the  Wiener  filter  for  wide-sense  stationary  and  un¬ 
correlated  object  and  noise.  This  exciting  result  is  ob¬ 
tained  without  knowledge  of  the  noise  psd  or  the  degra¬ 
dation  filter  H.  The  only  information  required  is  an  es¬ 
timate  of  the  object  psd  (Sff)  and  the  degraded  image 
itself.  However,  this  technique  yields  only  the  magni¬ 
tude  of  the  restoring  filter.  Unless  the  degrading  filter 
is  known  to  contain  no  zero  crossings  (e.g.,  Gaussian 
blur),  additional  methods  such  as  those  described  in  Sec¬ 
tion  3.2  must  be  used  to  estimate  the  restoration  filter 
phase.  Figure  7  shows  the  effectiveness  of  this  technique 
on  the  blurred  images  in  Fig.  3. 

4.2.6  Recursive  Filtering 

Recursive  filtering  is  a  specialized  restoration  tech¬ 
nique  that  assumes  that  the  original  object  is  corrupt¬ 
ed  only  by  additive  white  noise.  The  degradation  psf 
of  the  imaging  system  reduces  to  h(x,y)  =  6(x,y),  a 
two-dimemsional  dirac  function,  for  this  technique  and 

g(x,y)  =  f(x,y)  +  r>(x,y)  .  (61) 

The  object  f(x,y)  is  modeled  as  a  two-dimensional 
wide-sense  Markov  process.  Recall  that  for  a  one-dimen¬ 
sional  Mh  order  Markov  process,  the  value  of  the  func¬ 
tion  at  a  given  sample  location  is  dependent  only  on 
the  values  of  the  preceding  N  samples.  A  two-dimen¬ 
sional  (N,M) th  order  Markov  process  stipulates  that  the 
sample  value  at  a  given  location  is  dependent  on  the 
values  of  the  preceding  IVbyM  block  of  samples, 

Pr  { f(x  +  1,  y  +  1)  |  fUJ)  for  all  i  <  x,  j  <  y] 

=  Pr  i  f(x  +  1,  y  +  1)  |  fVJ)  for 

x  ~  N  <  i  <  x,  y  -  M  <  j  <  y\  .  (62) 


4*1 
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Figure  7  Homomorphic  deconvolution  restoration  of 
Fig.  3  with  (a)  SNR  =  33  dB  and  (b)  SNR  =  23  dB.’2 

In  image  restoration  applications,  this  technique  is 
frequently  used  assuming  that  the  original  object  is  a 
simple  two-dimensional  first  order  Markov  process.  The 
linear  estimate  of/(.v+l,  v+1)  has  been  found  to 


fix  +  l,v  +  I)  =  c/)f(x  -l-  l,r)  +  azf(x,y  -f  1) 

+  v)  +  «4  S  ( v,y )  .  (63) 


P*C(1,0)  +  p.C(O.l)  -  PxPyC(Q,0)  . 


C(0,0)  +  of 

where  C  is  the  object  and  estimate  cross  covariance, 
C(Ax,  Ay)  =  £([/U,y)  -  f(x,y)) 

[ fix  +  Ax,  y  +  Ay)  -  f(x  +  Ax,  y  +  Ay)]|  ,(65) 

and  a\  is  the  variance  of  the  white  noise. 

For  Eq.  64,  the  autocorrelation  function  of  the  ob¬ 
ject  is  postulated  to  be  spatially  separable  and  equal  to 

Rff  (Ax,  Ay)  =  p,  ^  pv  |iv|  .  (66) 


Recursive  filtering  is  substantially  different  from  the 
other  LSI  restoration  techniques  previously  described 
because  it  operates  in  the  spatial  and  not  the  Fourier 
domain.  Since  the  point  spread  matrix  of  the  filter  is 
of  limited  non-zero  extent  (four-element  array  for  a  first 
order  Markov  process),  the  spatial  convolution  can  be 
implemented  quickly  on  the  computer.  A  distinct  disad¬ 
vantage  of  this  approach  is  that  the  adjacent  sample 
correlations  must  be  assumed  constant  throughout  the 
entire  image  to  yield  a  spatially  invariant  restoration  fil¬ 
ter.  This  is  clearly  not  realistic  for  most  images,  which 
may  contain,  for  example,  regions  of  essentially  constant 
background  (px  =  py  =  1)  in  addition  to  regions  of 
rapidly  varying  intensities  (p*,pv  <<  1). 

This  gives  rise  to  the  development  of  the  spatially  vari¬ 
ant  restoration  technique  of  regional  recursive  filtering. 
Each  regional  filter  may  be  applied  to  areas  of  the  im¬ 
age  that  have  similar  spatial  correlation  properties.  Seg¬ 
mentation  of  the  image  into  these  regions  and  the  de¬ 
termination  of  the  individual  recursive  filters  clearly  add 
to  the  complexity  of  the  technique.  Regional  recursive 
filtering  also  creates  artifacts  at  the  region  boundaries 
that  must  be  smoothed  by  additional  processing. 
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4.3  LINEAR  SPATIALLY  VARIANT 
TECHNIQUES 

Section  4.2  outlined  those  restoration  techniques  ap¬ 
plicable  to  LSI  systems.  Because  the  imaging  system  was 
modeled  as  spatially  invariant,  the  powerful  theory  of 
Fourier  domain  analysis  could  be  used.  However,  real 
imaging  systems  are  seldom  spatially  invariant  and  are 
more  appropriately  modeled  as 

oc 

g  ( -v,  ,y, )  =  |  \  h(x,  ,  y,  ,x„  ,y„ ) 

-  * 

X  /(-Y„,v„  )  dx„  dy„ 

+  n(x,,\\)  ,  (67) 


coordinate  transformation  of  the  result  yields  the  re¬ 
stored  image  in  the  original  coordinate  system. 

The  application  of  this  approach  is  limited  to  those 
LSV  systems  that  are  decomposable  into  the  general 
form  shown  in  Fig.  8.  However,  the  effects  of  many 
forms  of  spatially  variant  image  degradation  can  be  de¬ 
scribed  by  this  decomposition.  The  applicability  of  this 
coordinate  distortion  technique  for  restoring  images  de¬ 
graded  by  certain  types  of  camera  motion 19,20  and  by 
the  optical  aberrations  of  astigmatism  and  curvature  of 
field*  has  also  been  shown. 

For  example,  a  simple  rotation  of  the  camera  about 
the  optical  axis  yields  the  LSV  result 

f  /(- W’») 

g(x,,y,)  =  \  — - —  ds„  (69) 

+  y"„ 


where  the  system  psf  is  a  function  of  both  the  object 
and  image  coordinates. 

The  simple  solutions  afforded  previously  by  Fourier 
analysis  techniques  no  longer  apply  because  of  the  spa¬ 
tial  variance  of  the  system  psf.  Restoration  techniques 
for  LSV  systems  are,  therefore,  fewer  in  number  and 
more  difficult  to  implement. 

4.3.1  Isoplanatic  Patches 

The  simplest  approach  to  L  SV  restoration  is  to  di¬ 
side  the  blurred  image  into  regions  or  isoplanatic 
patches.  In  each  isoplanatic  patch  the  spatially  invari¬ 
ant  assumption  is  approximately  valid  so  that  we  may 
model  the  system  as  piecewise  LSI, 


along  x2  +  y2  =  x?  +  y2  . 

The  limits  of  integration  are  from 

x„  =  .v,  x„  =  x,  cos  uiT  -I-  y,  sin  ef 
to 

y„  =  y,  y„  =  x,  sin  co7'  +  y,  cos  u>T 

Here, 

u)  =  constant  angular  rotational  velocity, 
T  =  time  of  photographic  exposure, 
ds„  =  differential  path  element, 


8,  (x,  ,y, )  =  \  \  h,  (x,  -  ,v„,  y,  -  y„  ) 

x 

x  ./(v„,y„)  dx„  dy„  ,  (68) 

for  (x,,\\  )  in  the  ,/rh  isoplanatic  patch.  Any  of  the 
previously  derived  spatially  invariant  restoration  tech¬ 
niques  can  then  be  applied  to  the  individual  patches. 
As  with  regional  recursive  filtering,  postprocessing  is 
usually  required  to  reduce  artifacts  generated  at  the  re¬ 
stored  image  isoplanatic  region  boundaries. 

4.3.2  Coordinate  Distortion  Method 

A  novel  approach  to  spatially  variant  image  restora¬ 
tion  developed  by  Sawchuk  and  Robbins  and  Huang6 
hinges  on  finding  a  nonlinear  coordinate  transformation 
that  maps  the  imaging  system  to  a  linear  spatially  in¬ 
variant  domain.  An  LSI  technique  can  then  be  used  to 
restore  the  image  in  this  new  domain,  and  an  inverse 


and  usT  is  constrained  to  be  less  than  2ir  radians  (over¬ 
all  camera  motion  must  be  less  than  one  full  rotation). 

The  imaging  system  is  otherwise  considered  to  be  per¬ 
fect,  and  no  noise  is  introduced  for  this  simple  exam¬ 
ple.  By  transforming  to  polar  coordinates  (r,  0),  we  ob¬ 
tain 


g(r,,e, ) 


1 

1 

I 


0 

-*■/ 


Ar,  ,.9„) 


CO 


d0„ 


f(r„A  -  0) 

GO 


c/0 


oc 

/( r„ ,  e,  -  0)  h(0)  dd  . 

--  oc. 


(70) 


'''A.  A.  Sawchuk,  “Space- Variant  Image  Motion  Degradation 
and  Restoration,"  Proc.  IEEE  60,  854-861  (1972). 

;"A.  A.  Sawchuk,  “Space- Variant  System  Analysis  of  Image 
Motion,"  J.  Opt.  Soc.  Am.  63.  1052-1063,  (19'’3). 
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Figure  8  LSV  coordinate  distortion  restoration. 


where  /;((()  =  (l/w)  for  -  u, >T  <  6  <  0.  The  system  is 
spatially  invariant  with  respect  to  the  new  coordinate 
domain.  An  LSI  technique  can  now  be  used  to  restore 
the  image  before  transforming  back  to  Cartesian  coor- 
dinates.  Figure  9  shows  the  usefulness  of  this  technique 
for  an  image  blurred  bv  camera  rotation.  The  LSI  tech- 


/  4 


\  j  • 


*  s 


nique  used  to  develop  the  restored  image  in  Fig.  9c  was 
a  simple  inverse  filter. 

The  drawbacks  of  this  technique  stem  from:  a)  being 
able  to  define  the  coordinate  transformation  for  a  given 
LSI  image  restoration  application;  and  b)  determining 
the  object  and  noise  stochastic  properties  in  the  new 
coordinate  system.  The  latter  point  is  relevant  when  LSI 
estimation  is  used  to  restore  the  degraded  image.  For 
example,  we  cannot  expect  that  minimizing  the  mean 
squared  estimate  error  in  the  transformed  domain  will 
necessarily  minimize  the  mean  squared  error  in  the  origi¬ 
nal  coordinate  system.  Flowever,  for  specific  applica 
lions  (e.g.,  camera  motion  and  aberration  correction), 
this  technique  may  be  quite  useful. 

4.3.3  Maximum  A  Posteriori  Restoration 

Maximum  a  posteriori  (MAP)  restoration  is  a  non¬ 
linear  iterative  technique  that,  in  its  most  general  form, 
is  applicable  to  spatially  variant  imaging  systems.  We 
describe  the  technique  for  a  discrete  imaging  system 
where  the  sampled  image,  object,  and  noise  data  are 
lexicographically  ordered  into  one-dimensional  vectors, 

g  =  IH\  f  +  n  ,  (71) 

and  His  a  two-dimensional  matrix  representing  the  ef¬ 
fects  of  the  spatially  variant  system  psf. 

The  MAP  estimate  of  the  vector  f  is  given  by 


Max  f/Mflg)) 

I  f  fv 


Figure  9  Coordinate  distortion  restoration  of  rotation¬ 
al  blur;  (a)  rotationally  blurred  object,  (b)  blurred  object 
transformed  to  pclar  cooidinates,  (c)  LSI  restoration 
of  transformed  object,  (d)  restored  image  transformed 
back  to  Cartesian  coordinates. 12 


Bayes's  rule  states 


Pr|f|gl  = 


Prlglfl  PrtJI 
Pr  Ig) 
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which  is  the  expression  we  need  to  maximize  with  re¬ 
spect  to  f .  If  the  original  object  and  noise  can  be  mod¬ 
eled  as  multivariate  Gaussian  processes,  i.e., 

/V{ f |  =  At  expf  -!(f-f)r  [Cyr1  (f  -  f)j  , 

L  2  )  (74) 

Pr\n\  =  A,  exp^-  ^nr  [C„] -1  nj  ,  (75) 

where  [Cf\  and  [C„  1  are  the  object  and  noise  covari¬ 
ance  matrices,  the  MAP  estimate  is  found  by  iteratively 
solving 

f\,AP  =  f  +  [Q]//r[Cn]-'  [g  -  //fMAp]  (76) 

for  fMAP .  Although  the  multivariate  Gaussian  assump¬ 
tion  may  seem  overly  restrictive  at  first,  it  gives  a  rea¬ 
sonable  model  for  many  images.  The  underlying  struc¬ 
ture  of  the  object  can  be  modeled  in  the  mean  vector 
f,  with  the  finer  structure  described  by  Gaussian  varia¬ 
tions  about  this  mean  vector.  This  may  be  a  more  realis¬ 
tic  model  than  that  obtained  by  assuming  a  wide-sense 
stationary  object.  Recall  that  the  stationarity  assump¬ 
tion  requires  that  every  element  in  the  mean  vector  f 
be  equal  to  the  same  constant  value  (as  in  Wiener,  ge¬ 
ometric  mean,  and  recursive  filtering). 

The  MAP  restoration  technique  requires  a  great  deal 
of  prior  knowledge  as  Eq.  76  demonstrates.  The  object 
and  noise  stochastic  properties,  as  well  as  the  imaging 
system  degradation  effects,  must  be  specified.  The  iter¬ 
ative  nature  of  this  technique  also  makes  it  computa¬ 
tionally  expensive  and  convergence  is  not  always  guaran¬ 
teed.  However,  MAP  restoration  has  been  used  success¬ 
fully  in  certain  applications21  and  can  be  extended  to 
nonlinear  spatially  variant  imaging  systems  as  well. 

4.3.4  Superresolution  Techniques 

The  limit  of  attainable  image  resolution  using  the  defi¬ 
nition  derived  from  the  Rayleigh  criterion2  is  deter¬ 
mined  by  the  cutoff  spatial  frequency  of  the  imaging 
system.  Even  with  the  best  possible  restoration  tech¬ 
nique,  object  information  outside  of  the  system  band¬ 
width  is  lost;  we  would  expect  a  minimum  resolvable 
spatial  element  corresponding  to  the  inverse  of  the  sys¬ 
tem  cutoff  frequency.  However,  if  certain  conditions 


!B.  R.  Hunt,  “Digital  Image  Processing,"  Proc.  IEEE  63, 
693-708  (1975). 


can  be  placed  on  the  object,  it  is  possible  (in  theory, 
at  least)  to  determine  the  object  spectrum  outside  of 
the  system  bandwidth.  The  restored  image  would  then 
achieve  the  property  of  spatial  superresolution — reso¬ 
lution  beyond  the  limit  imposed  by  the  system  cutoff 
frequency.  Two  techniques  that  have  this  property  are 
maximum  entropy  restoration4,5  and  the  method  of 
analytic  continuation.22,23  Unfortunately,  both  tech¬ 
niques  require  extremely  high  SNRs  to  achieve  super¬ 
resolution. 

Maximum  entropy  restoration  is  based  on  modeling 
the  original  object  vector  f  (suitably  normalized)  as  a 
probability  density  function  (pdf).  This  condition  guar¬ 
antees  positive  values  for  the  restored  image  vector  f. 
A  derivation  based  on  multinomial  pdf  (e.g.,  see  Refs. 
4  and  5)  yields 

f  Max  entropy  —  CXp  {  —  1  —  2 \[H  ) 

(g  —  //  f  Max  entropy  ))  ,  (77) 

where  the  additional  constraint  llg  -  //  f B 2  =  II n II 2 
has  been  incorporated  using  the  method  of  Lagrange 
multipliers.  Maximum  entropy,  like  MAP  restoration, 
is  a  nonlinear  iterative  technique.  The  estimate  derived 
is  intrinsically  not  band  limited  and  thus  superresolu¬ 
tion  results  in  high-SNR  cases. 

The  method  of  analytic  continuation  is  based  on  an 
eigenfunction  expansion  of  the  spatially  limited  object 
whose  Fourier  spectrum  is  known  only  in  the  system 
bandwidth.  Prolate  spheroidal  wave  functions  have  the 
interesting  properties  of  being  a)  orthonormal  on  the 
real  line  and  b)  a  complete  orthogonal  basis  for  the  set 
of  functions  band  limited  to  a  specified  bandwidth.  If 
the  imaging  system  is  modeled  as  a  spatially  separable 
system  defined  by  a  rectangular  aperture,  the  dual  char¬ 
acteristics  of  the  prolate  spheroidal  wave  functions 
(pswf)  lead  to  a  non-band-limited  estimate  for  the  origi¬ 
nal  object  (see  Refs.  22  and  23  for  details).  However, 
image  SNRs  on  the  order  of  30  dB  or  more  are  required 
for  analytic  continuation  to  result  in  useful  restoration. 
In  addition,  the  pswf  and  the  resulting  image  estimate 
are  difficult  to  compute,  which  further  limits  the  use¬ 
fulness  of  this  technique. 


22C.  K.  Rushforth  and  R.  W.  Harris,  “Restoration,  Resolu¬ 
tion,  and  Noise."  J.  Opt.  Soc.  Am.  58,  539-545  (1968). 
2,C.  L.  Rino,  “Bandlimited  Image  Restoration  by  1. inear 
Mean-Square  Estimation,”  J.  Opt.  Soc.  Am.  59,  547-553 
(1969). 
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5.0  COMPARISON  OF  RESTORATION  TECHNIQUES 


The  restoration  techniques  discussed  in  Sections  4.2 
and  4.3  are  summarized  in  Table  3.  As  the  second 
column  of  the  table  demonstrates,  the  LSI  restoration 
filters  appear  to  be  strikingly  similar.  We  have  seen  that 
these  filters  can  even  be  equivalent  under  certain  con¬ 
ditions;  homomorphic  deconvolution  is  equivalent  to 
geometric  mean  filtering  if  the  object  and  noise  are  un 
correlated  and  wide-sense  stationary.  Also,  constrained 
deconvolution  is  equivalent  to  parametric  Wiener  filter¬ 
ing  if  the  constraint  is  chosen  such  that  C  =  S„„/S ,j. 
Even  nonlinear  maximum  entropy  restoration  can  be 
related  to  a  simpler  LSI  technique.  If  a  linear  approxi¬ 
mation  to  the  exponential  in  Eq.  77  is  made  and  [//] 
is  assumed  LSI,  discrete  pseudoinverse  filtering  resto¬ 
ration  results. 

Given  the  similarities  between  the  restoration  filters, 
how  does  a  user  choose  one  technique  over  another  for 
a  particular  application?  There  are  important  factors 
that  differentiate  the  image  restoration  techniques  and 
act  as  guidelines  for  application,  including  the  following: 

1.  The  restored  image  quality  or  perceptibility  to  a 
human  observer  is  an  important  performance  fac¬ 
tor  that  depends  on  the  restoration  technique  used. 

2.  The  noise  performance  of  a  given  technique  is  a 
function  of  the  overall  image  SNR  as  well  as  the 
type  and  severity  of  the  image  degradation. 

3.  Each  technique  requires  certain  a  priori  knowl¬ 
edge  and  imposes  restrictive  assumptions  on  the 
imaging  system  model  which  limit  its  application. 

4.  The  feasibility  of  using  a  restoration  technique  for 
a  specific  application  is  often  determined  by  the 
length  of  time  required  to  implement  the  al¬ 
gorithm. 

Most  of  these  factors  are  delineated  in  Table  3.  One 
other  factor  that  influences  the  performance  of  the  resto¬ 
ration  technique  is  the  nature  of  the  system  noise.  The 
noise  type,  auto-  and  crosscorrelation  properties,  and 
action  on  the  system  (additive,  multiplicative)  will  also 
effect  the  performance.  This  effect  is  difficult  to  char- 
acteri  '  ■  other  than  by  stating  that  we  would,  of  course, 
expect  reduced  restoration  performance  for  techniques 
with  restrictive  system  assumptions  that  do  not  fit  the 
applicat,  ill. 

A  quick  (albeit  coarse)  comparison  of  the  most  com¬ 
monly  used  restoration  techniques  for  the  four  perfor¬ 
mance  criteria  defined  above  is  given  in  Table  4."1"4 

'T.  M.  Cannon,  H.  .1.  Trussell.  and  B.  R.  Hunt,  “Compar¬ 
ison  ol  Image  Restoration  Methods,"  Appl.  Opt.  17, 
33X4  3390  ( 1978). 


That  the  visual  quality  of  the  MAP,  homomorphic,  and 
geometric  mean  filters  is,  in  general,  better  than  Wien¬ 
er  filtering  may  be  surprising  at  first,  since  Wiener  filter¬ 
ing  supplies  the  MMSE  estimate.  However,  Wiener 
filtering  yields  the  MMSE  estimate  for  wide-sense  sta¬ 
tionary  object  and  noise  functions.  This  assumption  is 
seldom  true  for  real  images.  Also,  the  Wiener  filter 
trades  improved  noise  performance  for  reduced  image 
resolution.  As  studies  of  the  human  visual  system  have 
shown,25  human  observers  are  usually  willing  to  accept 
some  additional  image  noise  in  order  to  gain  improved 
resolution. 

The  noise  performance  of  the  techniques  shows  that 
inverse  filtering,  as  expected,  provides  unsatisfactory 
restoration  in  low-SNR  environments.  This  technique 
is  especially  sensitive  to  high-frequency  noise.  The  other 
restoration  filters  listed  in  Table  4  work  well  in  noisy 
environments.  As  the  SNR  increases,  however,  the  resto¬ 
ration  achieved  by  any  of  the  techniques  converges  to 
that  of  simple  inverse  filtering. 

The  degree  of  information  necessary  for  each  tech¬ 
nique  varies  widely  for  the  restoration  methods  listed. 
Homomorphic  deconvolution  requires  the  least  amount 
of  a  priori  knowledge  while  Wiener,  geometric  mean, 
and  MAP  filtering  require  detailed  information  about 
the  imaging  degradation  and  stochastic  properties  of 
the  object  and  noise.  The  degree  of  information  need¬ 
ed  for  inverse  filtering  and  constrained  deconvolution 
falls  between  these  two  extremes. 

Often,  the  most  important  factor  in  choosing  a  resto¬ 
ration  technique  is  its  computational  complexity.  If  large 
amounts  of  image  data  need  to  be  processed  quickly 
and  cost  effectively,  it  would  be  implausible  to  use  an 
iterative  algorithm  (without  guaranteed  convergence)  like 
MAP  restoration.  Even  :he  partitioning  required  to 
compute  and  apply  the  homomorphic  filter  may  take 
too  much  time  for  such  an  application.  Inverse,  Wiener, 
geometric  mean,  and  constrained  deconvolution  filter¬ 
ing  have  moderately  fast  implementations  using  two- 
dimensional  fast  Fourier  transform  (FFT)  algorithms. 
Wiener  filtering  has  an  additional  advantage;  the  resto¬ 
ration  filtering  can  be  performed  in  any  unitary  trans¬ 
form  domain.  Fourier,  Hadamard,  identity,  or  Karhun- 
en-Loeve  transform  processing  is  chosen  depending  on 
which  technique  yields  the  fastest  implementation  for 


'  1 .  G.  Stockham,  .Ir . ,  “Image  Processing  in  the  Context  of 
a  Visual  Model,"  Proc.  ,'LLf:  60,  828-842  (1972). 
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Table  3 

Summary  of  image  restoration  techniques. 


Restoration  Restoring 

technique  filter 


Restoration 

criterion 


A  priori  knowledge 
and  assumptions 


Performance 

advantages 


Disadvantages 


Restoring  filter 
variations 


Inverse  filtering  .Vf  = 


Minis  -  h*/): 


LSI  system 
H  known 


Computational 

simplicity 


Amplifies  high- 
frequency  noise 
Indeterminate  when 
H  =  0 


1  [  -S„„ 

Wiener  tillering  A f  -  -  II  +  — r— 

HI  \HS„ 


Min  £||/  -  /| :  1  LSI  system 

H.  S„„,  S',  known 
n,  /  stationary 


1  S„„  V 

Geometric  mean  A/  =  I  +  7  — 

_ _  H  '  W|  S.,  J 


filtering 


LS!  system 
H,  S,„.  S„  known 
n  /  stationary 


Combines  good  low- 


I  f 

C  onstramcd  \t  =  l  *  . 

.  I,».  \nv  nil  it  mn  ^ 


Min  |r*/;: 
subject  to 

I*  *  •  A* 

=  o; 


LSI  system 
H,  o;  known 
C  user  defined 


User -defined  constraint 
generates  desired 
restoration  without 
knowledge  of  S/f,  S„„ 


Homomorphic  V/ 
deconvolution 


2>' 


s„  =  -V 


LSI  system 
Sfl  known 

Extent  of  H  <  area  of 
partition  block 
Ergodic  g 


Recursive 
lilt  enng 


Mm  £li/  -  /I  *  |  Ideal  I.SI  imaging  system  Quick  filter 

White  additive  noise  implementation 

f{x.v)  is  (\.Af)th  order 
wide-sense  Markov 
R„  and  C„  known 


Iw>planafic  parch  Any  LSI  filter 
filtering 


I.SV  system  is  piecewise 
LSI  plus  LSI  technique 
assumptions 


Straightforward 

implementation 


<  oordmate  Any  LSI  Hirer 

distortion  merhixl 


LSV  system  described  as 
I.SI  after  known 
coordinate  transforma¬ 
tion  plus  LSI  technique 
assumptions 


Fast  LSV  restoration 
achieved  when 
coordinate  transforma¬ 
tion  known 


MAP'  restoration 


Max  Pr[  fig) 


LSV  system 
/  and  f>  multivariate 
Gaussian  processes 
[«l.  L  IG|.  |C„| 
known 


Visually  pleasing 
restoration 


Maximum  ‘ 
enirops 


f  In  fj  L.SV  system 

f,  f  non-negative 
| //].  n'„  known 


Visually  pleasing 
restoration 
Supcrresolution 


AnaMic 

continuation 


Mm  fl./  t:  I 


1SV  system 
/  analytic  and  band- 
limitcd  to  square  band 
region 

Prolate  spheroidal  wave 
functions  known 


Supcrresolution 


for  details  regarding  ihis  technique,  see  Section  4  2  ft 
See  Section  4  \  \ 

See  Section  4  t  4 
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Pseudoinverse 
Af  =  H*  m 


,'Lmo[yrw] 


Optimum  linear  MSE 
estimate 

Good  noise  performance 
M  =  0  when  H  =  0 


Stationary  /  assumption  Parametric  Wiener 


unrealistic 
Restored  image  is 
smoothed 


M  = 


frequency  restoration 
of  inverse  filtering  with 
high-frequency  noise 
performance  of  Wiener 
filtering 

User  tailored  for  each 
image 


Stationary  /  assumption  True  geometric  mean, 


unrealistic 

Iterative  and  subjective 
method  to  find  param¬ 
eters  y,  a 


a  =  lA,  y  =  I 
Inverse  filter,  a  -  I 
Parametric  Wiener, 
a  =  0 


Iterative  method  to 
find  y 

Restoration  more  noise 
sensitive  than  Wiener 
filtering 


Finite  second  degree 
difference  matrix 
|c(ij)| 

Human  visual  response 
model  (c(i'j)| 


No  knowledge  of  H 
required 

Good  noise  performance 


Magnitude-only  restora¬ 
tion  filter  requires 
additional  processing  to 
determine  filter  phase 
Computationally  complex 


Equivalent  to  geometric 
mean  with  a  = 

7  =  I  for  stationary  /, 


Restrictive  and  unrealistic  LSV  regional  recursive 
assumptions  filtering 


Computationally 

expensive 


Only  works  for  certain 
decomposable  LSV 
systems  (Fig.  8) 


Iterative  technique 
Convergence  not 
guaranteed 


Iterative  technique 
Convergence  not 
guaranteed 

Very  susceptible  to  noise 


Computationally 

expensive 

Requires  very  high  SNR 
( >  30  dB) 
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Table  4 

Comparison  of  restoration  techniques. 

Visual  quality  Noise  Degree  of  required  Computational 

(moderate  SNR)  performance  a  priori  knowledge  complexity 


Inverse  filtering 

Fair 

Poor 

Low 

Low 

Wiener  filtering 

Good 

Good 

High 

Moderate 

Geometric  mean 
filtering 

Very  good 

Good 

High 

Moderate 

Constrained 

deconvolution 

Good 

Fair 

Moderate 

Moderate 

Homomorphic 

deconvolution 

Very  good 

Good 

Low 

High 

MAP  restoration 

Very  good 

Good 

High 

Very  high 

the  particular  image.26  (See  Ref.  12  for  restoration  al¬ 
gorithms  and  additional  insight  on  their  implementa¬ 
tion.) 

Ultimately,  the  choice  of  restoration  technique  is  de¬ 
termined  by  the  particulars  of  the  application.  What 
kind  of  visual  quality  is  required?  A  restored  image  that 
will  be  subsequently  processed  by  additional  software 
for  classification/detection  purposes  will  not  necessar¬ 
ily  require  the  same  degree  of  visual  quality  as  one  that 
will  be  presented  to  a  human  observer.  What  is  the  sys¬ 
tem  SNR  and  how  severe  is  the  degradation?  If  the  sys¬ 
tem  consistently  operates  at  high  SNRs,  simple  inverse 
filtering  may  suffice.  What  kind  of  prior  knowledge  is 
available?  If  little  information  about  the  system  is 
known  or  estimable,  homomorphic  deconvolution  may 
be  the  most  attractive  technique.  However,  the  user  must 
remember  the  assumptions  couched  in  this  method:  the 
non-zero  extent  of  H  is  less  than  the  partition  size  and 
H  is  real  and  non-negative. 


And  finally,  what  are  the  computational  considera¬ 
tions  for  the  particular  application''  Are  there  reels  of 
data  or  a  single  image  to  restore?  Is  the  data  stored  and 
processed  in  an  image  processing  laboratory  or  must 
restoration  be  accomplished  in  close  to  real  time?  For 
very  fast  restoration  applications,  the  format  of  the  im¬ 
age  data  can  further  influence  the  choice  of  restoring 
technique.  For  example,  video  image  data  received  one 
scan  line  at  a  time  lends  itself  to  techniques  that  can 
be  implemented  by  one-dimensional  line  processing. 
Even  though  faster  two-dimensional  FFT  algorithms  ex¬ 
ist  that  operate  on  a  two-dimensional  extension  of  the 
one-dimensional  Cooley-Tukey  butterfly  operation,27 
most  two-dimensional  FFT  algorithms  are  performed 
by  sequential  one-dimensional  FFT  computations.  This 
makes  techniques  such  as  inverse,  Weiner,  geometric 
mean,  and  constrained  deconvolution  filtering  especially 
applicable  to  scanned  imagery. 


>VV.  K.  Pratt,  “Generalized  Wiener  Filtering  Computation 
Techniques,”  IEEE  Tram.  Comp.  C-21,  636-641  (1972). 


27D.  B.  Harris,  J.  H.  McClellan,  D.  S.  K.  Chan,  and  H.  W. 
Schuessler,  “Vector  Radix  Fast  Fourier  Transform,”  IEEE 
Conf.  on  Acoustics,  Speech,  and  Signal  Processing,  Hart¬ 
ford,  Conn.  (May  9-11,  1977). 
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6.0  CONCLUSIONS 


Although  the  techniques  discussed  above  represent 
those  most  frequently  used,  they  are  only  a  subset  of 
the  image  restoration  techniques  that  have  been  devel¬ 
oped.  Some  additional  methods  are  maximum  likeli¬ 
hood,  maximum  weighted  Burg  entropy,  and  even  mini¬ 
mum  entropy.12  Frieden28  unifies  Bayesian  estimate 
methods  under  a  general  theory  of  maximum  physical 
likelihood  (note,  however,  that  Frieden’s  definition  of 
maximum  likelihood  differs  from  convention).  Several 
nonlinear  and  linear  recursive  restoration  techniques  are 
studied  in  another  publication  by  Meinel.29 

With  the  exception  of  maximum  entropy  and  MAP 
restoration,  the  techniques  in  Section  4.0  do  not  guar¬ 
antee  non-negative  values  for  the  restored  image,/.  If 
the  optical  intensity  of  the  object  is  measured  by  the 
function  /  (i.e.,  /  >  0),  we  would  prefer  that  the  re¬ 
stored  image  be  non-negative  also  (f  >  0).  This  can  be 
done  by  adding  a  positivity  constraint  to  any  restora¬ 


tion  criterion.  Unfortunately,  this  usually  does  not  lead 
to  a  closed  form  expression  for  the  restored  image  and 
requires  computationally  intensive  nonlinear  program¬ 
ming  to  obtain  a  solution. 

New  image  restoration  techniques  can  be  developed 
simply  by  defining  new  restoration  criteria.  Most  current 
techniques  are  based  on  criteria  that  generate  visually 
pleasing  restorations.  However,  if  the  restored  image 
is  subsequently  processed  by  additional  software  rather 
than  presented  to  a  human  observer,  it  may  be  possible 
to  develop  new  techniques  that  optimize  the  perfor¬ 
mance  of  this  postrestoration  processing.  One  applica¬ 
tion  where  this  may  prove  profitable  is  in  the  area  of 
automatic  scene  detection  and  classification.  Potentially, 
one  could  develop  a  restoration  technique  based  on  a 
criterion  that  would  be  tailored  to  give  maximum  prob¬ 
ability  of  detection  and  correct  classification  of  the  resul¬ 
tant  restored  image  scenes. 


'B.  R.  Fricden,  “Unified  Theory  for  Estimating  Frequency- 
of-Oceurrence  Laws  and  Optical  Objects,”  J.  Opt.  Soc.  Am. 
73,  927-938,  (1983). 

T.  S.  Meinel,  “Origins  of  Linear  and  Nonlinear  Recursive 
Restoration  Algorithms,”  J.  Opt.  Soc.  Am.  3,  787-799 
( 1986). 
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